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i.  introduction 

This  paper  is  concerned  with  the  statistical  assessment  of  dependency 
beyond  autocorrelation  in  the  context  of  nonlinear  time  series  models.  The 
central  theme  is  that  residuals  from  models  fitted  according  to 
inapplicable  linearity  assumptions  can  also  profitably  be  used  for  further 
analysis  beyond  linearity.  A  well  known  property  of  linear  time  series 
residuals  is  that  they  are  uncor related;  in  the  context  of  linear  models 
they  should  also  be  independent,  apart  from  the  effects  of  parameter 
estimation  with  short  series.  However,  in  the  context  of  nonlinear  models 
it  is  not  often  recognised  that  uncorrelated  residuals  also  hold 
information  concerning  higher  order  dependence  in  the  data.  Developments 
of  such  a  higher  order  residual  analysis  will  be  explored  here,  in 
particular  for  two  types  of  nonlinear  autoregressive  model  which  have  the 
usual  linear  Yule-Walker  autoregressive  correlation  structure.  Some  higher 
order  dependency  correlations  will  be  obtained.  In  a  companion  paper,  the 
suggestion  of  reversed  residuals  will  be  made,  and  the  analysis  given  in 
this  paper  will  be  extended  to  encompass  these  reversed  residuals. 

2.   AUTOREGRESSION  AND  LINEAR  AUTOREGRESSIVE  RESIDUALS 
2.1  Auto regress ion  for  Nonlinear  Models 

The  standard  form  of  autoregression  needs  widening  for  use  with 
nonlinear  models;   we  consider  first  the  standard  form  which  is  explicitly 
autoregressive  in  a  linear  additive  way  and  then  several  weaker  variants. 
A  stationary  time  series  {Xt}  of  mean  ju  is  assumed.   Under  the  linear 
autoregressive  model,    of  order  p,  the  {Xt)  satisfy  the  equation 

Xt-^a1(Xt_1-^)+a2(Xt_2-/x)+...+ap(Xt_p-M)+et,    t=0,  ±1,  ±2,  .  .  .  (2.1) 
where  the  et  are  independent  and  identically  distributed  and  n,a1#a2, . . . ,ap 
are  fixed   parameters.   A  more  general  definition  of  autoregression  of  order 
p,  could  be  the  linear  conditional  expectation   requirement  that 
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E(Xt-ll\Xt-ltXt-2 Xt_p) 

=  a1(Xt_1-M)+a2(Xt_2-M)+...+ap(Xt-p-M)/    t=0,±l,±2 (2.2) 

The  definition  (2.1)   implies  (2.2)  but  not  vice-versa.    Thus  this 

definition  could  apply  to  models  which  are  not  of  the  linear  form  (2.1), 

either    because    the    e^  are    dependent,    but    still    with 

E(  €-t  |Xt-i,JCt_2/  .  . .  ,Xt_p)=0,  or  because  the  model  has  some  other  structure 

altogether.    For  instance,  there  are  the  random  coefficient  models  of 

Nicholls  and  Quinn  (1982),  the  exponential  distribution  random  coefficient 

models  of  Lawrance  and  Lewis  (1981,1984),  the  discrete  distribution  random 

coefficient  models  of  Jacobs  and  Lewis(1983),  and  the  gamma-beta  random 

coeffxcient  models  of  Lewis  (1981). 

Random  coefficient  autoregressive  models  of  order  p  take  the  general 

form 

Xt  =  At<1>  Xt_1+At(2)xt_2+...+At<P>Xt_p+Bt,     t=0,±l,±2 (2.3) 

where  the  vector  of  coefficients  {A-j-( 1  K&t^  2  ),  .  .  .  ,At-(P),B-t}  is  a  stationary 

vector  sequence  of  independent  random  variables,  and  sometimes  in  addition, 

Et  is  independent  of  the  random  coefficients  at  time  t.   It  is  easy  to  see 

that  such  models  satisfy  the  linear  conditional  expectation  definition  of 

pth  order  autoregression,  but  are  nonlinear.   A  first  order  case  of  the 

type  (2.3)  will  be  used  to  illustrate  the  proposed  method  of  residual 

analysis  of  autoregressive  nonlinearity. 

A  further  and  weaker  definition  of  autoregression  is  the  requirement 

that  the  autocovariances  of  the  {X^},  denoted  by  {y^},  just  satisfy 

Yule-Walker  linear  difference  equations  of  the  form, 

Vr=aiyr-l+a2yr-2+ • • • +apyr-p»  yr=V-r >  r^1 » 2 >  • • •  (2.4) 
for  suitable  constants  ai,<X2,  . . .  ,otp;  this  will  be  referred  to  as 
Yule -Walker  autor egression.  It  is  true  for  processes  which  satisfy  (2.1) 
and  (2.2),  as  may  be  verified  in  the  usual  manner  by  multiplying  X^  by  X^-r 
and  taking  expectations.  The  reverse  is  not  true,  (2.4)  does  not  imply 
(2.1)  or  (2.2).  A  case  in  which  (2.4)  holds  but  in  which  (2.1)  and  (2,2) 
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do  not,  is  the  product  autoregression  model  of  McKenzie  (1982),  where,  with 
p=l,  E(XtlXt-i)  is  a  fractional  power  of  x^-i;  this  model  will  also  be  used 
in  Section  4.3  as  an  ilustration  of  the  proposed  residual  methods. 

A  variety  of  other  models  can  satisfy  the  Yule-Walker  definition  of 
autoregression,  and  yet  not  satisfy  the  linear  conditional  expectation 
definition.  Amongst  these  are  first  order  Markov  chains  under  fairly 
general  conditons,  and  others  such  as  the  semi-Markov  generated  processes 
discussed  in  Cox  and  Lewis  (Chapter  7,  1966)  and  Lewis  (1980).  Notice  that 
we  do  not  define  nonlinear  autoregressive  models  in  a  constructive  way. 
The  class  is  so  wide  as  to  make  this  impossible;  one  such  class  has  been 
studied  by  Jones  (1978).  Rather,  we  require  that  the  autocorrelations 
should  satisfy  linear  equations,  similar  in  structure  to  those  satisfied  by 
the  autocorrelation  of  linear  autoregressive  models.  In  view  of  this,  our 
suggested  analysis  extends,  rather  than  superceeds,  conventional  methods. 

2.2  Definition  and  Discussion  of  Linear  Autoregressive  Residuals 

For  the  analysis  of  time  series  data  involving  models  satisfying  (2.1), 
or  (2.2)  or  minimally  (2.3),  the  use  of  linear  autoregressive  residuals  of 
order  p,  defined  as 

Rt<  P  >  -  (  Xt-u  )-<xx(  Xt-irll  )-<x2(  Xt-2-M  )"  •  •  •  -«p(  Xt-p-M )         (2.5) 
is  suggested.   This  suggestion  is  based  on  the  following  theorem,  which  is 
a  generalization  of  a  result  given  in  Lawrance  and  Lewis  ( 1984,  Section 
7.2)  for  p=2. 

Theorem.   Let  the  stationary  process  (Xt>  satisfy  the  Yule-Walker  type 
equations  (2.4).   Then  the  linear  autoregressive  residuals  {Rt*p)}  defined 
at  (2.5)  are  uncorrelated  (although  not  necessarily  independent). 
Proof.   The  autocovariances  of  the  residuals  (2.5)  are 
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Cov[Rt<  P  > ,  Rt+r(  P  }  ]  =  Cov[  (  Xt-ji ) ,  Rt+r(  P  >  ]-a!Cov[  (  Xt-i-M  ) ,  Rt+r(  P }  ]"  •  •  • 

-apCovtCXt-p-M-J^Rt+^P^  (2.6) 

=  Cov[(Xt-/i),Rt+r(P)]-a1Cov[(Xt-M),Rt+r+l(P)]--  •  • 

-apCov[(Xt-M),Rt+r+p(P)]-  (2.7) 

Equation  (2.7)  follows  because  the  {X^}  process  is  stationary  and 
consequently  the  {R^P^}  process  is  stationary.  The  covariances  in  (2.7) 
need  only  be  considered  for  positive  lag  since  the  autocovariance  is  an 
even  function  of  r.  Then  the  crosscovariances  on  the  right-hand  side  of 
(2.7)  are  all  of  the  same  type  and  given  by 
Cov[(Xt-/i)/Rt+r(P)]  =  Cov[(xt-M),{(xt+r-M )-<*!( xt+r_i-M).  .  . -<xp( Xt+r-p-/x ) }  ] 

=  yr_aiyr-l-- • --apYr-p      r=l, 2, ....         (2.8) 
Now  by  the  Yule-Walker  equations  (2.4),  the  expression  (2.8)  is  zero.   Thus 

using  (2.8)  in  (2.7) 

Cov[Rt(P>,Rt+r(P)]  =  0       r=±l,±2,...,      (2.9) 
as  was  to  be  proven.   The  proof  is  immediate  for  the  linear  autoregressive 

model  (2.1),  since  Rt^P^  =  «f  an<^  €t  an<*  Rt+r^P^  are  DY  definition 
independent . 

Note  that  the  linear  autoregressive  residuals  {Rt^P)}  will  generally 
still  be  dependent  though  uncorrelated  in  nonlinear  modelling  of  the  {X^}; 
with  the  linear  autoregressive  model  (2.1)  the  {Rt^P^}  will  not  only  be 
uncorrelated  but  also  independent.  It  is  this  difference  which  uill  be 
exploited  to  explore  nonlinearity  in  pth  order  nonlinear  autoregressive 
processes.  The  dependency  attributable  to  parameter  estimation  is  taken  to 
be  small  in  the  large  scale  applications  we  have  in  mind;  indeed,  nonlinear 
modelling  of  short  series  of  data  may  well  be  hard  to  justify. 

The  quantities  R^P)  are  autoregressive  residuals  in  the  sense  that 
r^(P)  is  the  residual  of  X^  after  subtracting  off  ji^P),  its  best  linear 
least  squares  predictor  in  terms  of  Xt-i,X^-_2,  .  .  .  ,Xt-p,  given  by 

jtt<  P  >  =  ^+ax(  Xt-i-ji  )+a2(  Xt_2-M  )+ .  .  .  +ap(  Xt-p-ji ) .  (  2  .  10  ) 

Thus  the  residuals  (R^P^)  give  the  basic  way  of  taking  out  the  linear 
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correlation  component  in  models  with  pth  order  autoregression  minimally  of 
the  'Yule-Walker'  form  (2.4).  Note  too  that  if  the  process  is 
autoregressive  in  this  sense,  linear  autoregression  on  terms  beyond  X^-p-M/ 
for  example  (Xt_p_x-/i),  will  yield  ap+k=0,  k=l,2,...;  this  is  a  property  of 
the  linear  component  being  taken  out  so  as  to  minimize  the  expected  mean 
square  of  the  residual  Rt(P+K). 

A  further  point  worth  noting  about  the  residuals  Rt/P^  concerns  their 
crosscovariances  with  the  X^'s.   In  the  proof  of  the  theorem  it  is  seen 
that  Cov[Xt-H,Rt+r^P*l  =  °  for  r  =  1,2,...  .   However,  the  other  half  of 
these  crosscovariances  is  non-zero. 
Example:  The  nonGaussian  linear  AR(  l )  model 

Taking  (2.1)  with  p=l,and  p  instead  of  a^,  the  crosscorrelation  of  X^-^  and 
Rt-r(P)  in  the  AR(1)  model  is  given  by 

c:orr(xt-/x,Rt_r)  =  (1-p2)1/2  pr  for  r=l,2,3,...;  (2.11) 
the  superscript  has  been  dropped  from  Rt-r*  as  *■*  wiH  b6  i-n  similar  future 
use.  Note  further,  that  in  this  case  Rt+r  is  independent  of  X^  for 
r=l,2,...  .  Further  use  of  (2.11)  will  be  made  in  Sections  3  and  5. 

For  the  use  of  {Rt^pM  in  data  analysis,  the  order  p  of  the  linear 
aspect  of  the  autoregression  needs  to  have  been  chosen;  any  of  the 
available  standard  methods  may  still  be  used.  In  addition,  of  course,  the 
jt,ai,a2, . .  .  ,<Xp  need  to  be  estimated;  there  are  at  least  two  convenient 
possibilities:  (1)  the  assumption  of  a  linear  autoregressive  model  like 
(2.1)  and  the  use  of  least  squares  estimation,  and  (2)  a  non-model  based 
approach  to  estimation  employing  the  first  p  Yule-Walker  type  equations 
(2.4).  The  latter  is  suggested  here;  however,  Tjostheim  and  Paulsen  (1983) 
recommend  (1)  when  dealing  with  modest  sized  samples  from  linear  models,  on 
account  of  serious  estimation  bias  with  the  Yule-Walker  estimates.  Since 
the  use  envisaged  here  is  primarily  with  nonlinear  models  and  substantial 
sets  of  data,  the  superiority  of  ( 1 )  over  ( 2 )  is  not  established  or 
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crucial . 

3.   ASSESSMENT  OF  HIGHER  ORDER  DEPENDENCE 

Since  the  linear  autoregressive  residuals  {R^P^}  of  (2.5)  have  zero 
covariances  when  (2.4)  holds,  but  for  nonlinear  processes  need  not  be 
independent,  a  residual  analysis  of  this  nonlinear ity  can  be  based  on  an 
assessment  of  their  higher  order  dependence.  To  consider  what  form  this 
might  take,  we  note  that  the  use  of  {X^2}  has  been  suggested  by  Granger  and 
Andersen  (1978, p. 63)  for  bilinear  models,  for  which  many  of  the  simple 
cases  have  zero  autocorrelations  themselves.  The  corresponding  initial 
suggestion  here  is,  as  was  briefly  illustrated  in  Lawrance  and  Lewis 
(1984),  to  use  the  residuals  {(Rt^^)2}-  Displays  can  then  easily  be  made 
of  the  associated  autocorrelation  functions,  scatter  plots,  periodograms, 
cumulative  periodograms,  etc,  using  standard  (second  order)  time  series 
software . 

However,  the  autocorrelations  of  {(Rt^^)2}  are  fourth  order  quantities 
in  the  original  series  {X^},  which  is  a  double  jump  from  the  second  order 
autocorrelations  of  the  series  {X-j-}.  Such  quantities  will  be  very 
difficult  to  handle  theoretically  with  most  types  of  nonlinear  model. 
The  crosscorrelation  function  of  {Rt^^}  and  {(Rt^^2)  is  essentially  3rd 
order,  but  by  involving  the  variance  of  {(Rt^^)2}*  also  needs  some  fourth 
order  join  I  moments  of  {X^-}  up  to  lag  p.  The  necessary  calculations  of  all 
these  quantities  for  two  nonlinear  models  will  be  given  in  Section  4. 

The  behaviour  of  the  crosscorrelation  function  of  R^P)  and  {(R^P^)2} 
may  be  judged  against  the  fact  that  the  {R^P^}  are  Independent  for  the 
linear  autoregressive  model  of  order  p,  and  hence  it  will  be  zero  except  at 
lag  zero.  For  the  random  coefficient  autoregressive  processes  considered 
in  Section  4,  the  suggested  crosscorrelation  function  will  be  shown  to 
posess  a  useful  cut-off  property;  this  generalizes  the  NEAR( 2 )  result  given 
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in  Lawrance  and  Lewis  (1984). 

Two  rather  simpler  crosscorrelation  quantities  can  be  proposed  from 
consideration  of  the  zero  covariance  result  of  (2.8),  and  the  additional 
fact  that  when  the  pth  order  linear  autoregressive  model  holds,  Xt  and 
Rt+r(P)  (r=l,2,...)  will  be  independent,  and  not  just  uncor related. 
Working  in  terms  of  {X^}  adjusted  for  its  mean  \i,  which  is  better 
computationally  and  often  nicer  theoretically,  the  following  may  be 
considered 

Corr[(Xt-/02,Rt+r(P)]'     for  r=0,±l, ±2, . . .  (3.1) 

Corr[Xt-/x,(Rt+r(p))2]'    for  r=0±l, ±2, . . .  .  (3.2) 

These  autocorrelations  are  not  equivalent  in  the  aspects  of  higher  order 
dependency  of  tho  {X^}  process  which  they  assess.  To  see  this,  note  that 
the  covariance  corresponding  to  (3.1)  involvon  only  third  order  joint 
moments  of  the  form  E[(Xt-*i)2(Xt+r-i-M)]/  i=0,+l, . . . , ±p  while  the 
covariance  corresponding  to  (3.2)  involves  additional  joint  moments  of  the 
form  E[Xt-^.)(Xt+r~^)(xt+r-i)]'  i=0, ±1, . . . , ±p.  Further,  the  denominator  of 
(3.2),  by  needing  var{(R^^P^ )2}/  is  more  complicated  in  the  higher  order 
moments  it  involves  relative  to  the  denominator  of  (3.1)  which  needs 
var{Rt(P)}.  However,  both  correlations,  and  particularly  (3.1),  are  more 
tractable  than  the  autocorrelations  of  {(Rt^p))2}  or  the  crosscorrelations 
of  {Rt<P>}  and  {(Rt(f>>)2}. 
Example  (Continued  from  (2.11)):   The  nondaunnian  linear  AR(1)  model. 

As  an  illustration  of  the  use  of  (3.1)  and  (3.2),  for  the  linear  AR(1) 

model,  ((2.1)  with  p=l  and  p  instead  of  a^ )  there  are  the  results 

0  r=-l,-2,... 


Corr[(Xt->i),(Rt_r)2]  = 


_-n  \2 


OorrCCXt-fO^Rt-r]    = 


skew(  X )  — ^ii ,( 1-P3  )pr ,     r=0 , 1 ,  .  . . 
var(R2) 

r=-l,-2, . . . 

(3.4) 


lskew(X)[Var{(x_^)2}]i7T  ^J   P  .  r-0,1,... 
Note  the  faster  geometric  decrease  in  (3.4)  relative  to  (3.3);  (3.4) 

contains  more  higher  moment  information.  The  other  halves  of  both  these 
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crosscor relation  functions  are  of  course  zero.  Behaviour  of  this  type  for 
these  functions  -  zero  for  r<0  and  geometric  decay  for  r>0  -  would  suggest 
linear  models  with  skewed  marginal  distributions  or  nonlinear  models. 
Discussion  of  parallel  results  for  two  types  of  nonliruMr  modol  are  given 
in  section  4. 

The  urn*  of  squaring  in  the  construction  of  these  higher  order 
dependency  measures  is  recognized  as  being  pragmatic  and  somewhat  arbitary; 
it  does  however  lead  to  expressions  involving  selected  types  of  simpler 
higher  order  cross  moments.  The  end  use  of  the  higher  order  dependency 
measures  can  either  he  exploratory  for  a  given  data  set,  to  ascertain 
Whether  there  is  appreciable  nonlinearity  present,  or  constructively,  to 
provide  evidence  for  fitting  suitable  types  of  nonlinear  model  which  can 
match  the  observed  higher  order  dependency. 

Earlier  discussion  of  higher  order  dependence  in  nonc'aussian  linear 
models  is  given  by  Rosenblatt  ( m80 )  in  terms  of  the  bispectrum;  aspects  of 
nonlinearity  and  higher  order  spectra  are  briefly  considered  in  Rosenblatt 
(1979). 
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4.   HIGHER  ORDER  DEPENDENCY  FOR  TWO  TYPES  OF  NONLINEAR  MODEL 

In  this  section  we  obtain  properties  of  the  measures  of  higher  order 
dependency  discussed  in  the  previous  section,  for  two  specific  types  of 
nonlinear  models.    The  models  considered  are  autoregressive  in  the 
Yule-Walker  sense  of  (2.4),  but  not  in  the  linear  sense  of  (2.1). 
4.1  Random  Coefficient  Autoregressive  Models 

As  already  remarked  in  Section  1,  a  generalization  of  the  linear 
autoregressive  model  (2.1)  is  to  let  the  coefficients  ai,<X2,...ap  be  random 
variables.  One  general  class  of  such  models  has  been  discussed  by  Nicholls 
and  Quinn  (1982)  who  cite  Andell  (1976)  and  articles  in  the  economic 
literature.  Other  classes  of  models  with  random  coefficients  include  the 
discrete  distribution  models  of  Jacobs  and  Lewis  (1983)  and  the  exponential 
models  of  Lawrance  and  Lewis  (1981,1984).  The  class  of  random  coefficient 
autoregressive  processes  to  be  of  interest  here  is  given  by 

Xt  «  At<1)xt_1+At(2)xt_2+...+At(P)xt_p+Bt        t=0,±l,±2 ,    (4.1) 

where  {A^1  ),Afc( 2  K  . .  •  ,Afc(P>,B.fc}  is  a  stationary  vector  sequence  of 
independent  random  variables  with  E(At(3))=<Xj  for  j=l,  .  .  .  ,p.  The 
components  of  the  vectors  are  not  necessarily  independent.  For  example  the 
discrete  distribution  models  of  Jacobs  and  Lewis  (1983)  can  be  written  in 
this  form,  and  have  dependent  coefficients,  as  do  the  exponential  models  of 
Lawrance  and  Lewis  (1981,  1984).  It  is  easily  verified  that  the  process 
(4.1)  satisfies  conditional  expectation  autoregression  (2.2)  and  thus  also 
the  weaker  Yule-Walker  definition  (2.4);  it  will  also  clearly  be 
stationary.  Note  also  that  the  standard  linear  autoregressive  model  (2.1) 
is  a  special  case  of  (4.1)  in  which  the  random  coefficients  are  actually 
constant. 

We  now  give  a  characteristic  result  for  this  type  of  process  when 
higher  order  dependency  is  measured  by  Corr[Rt<P),(Xt4.r-M)2]  or 
Corr[Rt<P>,(Rt4.r(P>)2],  assuming  that  the  Rt(P>'s  are  uncorrected . 
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Theorem.   With  the  random  coefficient  model  (4.1), 

Corr[Rt<P>,(Xt_r-ji)2]  and  Corr[Rt(P>,(Rt_r<P)  )2]  are  equal  to  zero  for 

Proof:   Using  the  definitions  (2.5)  and  (4.1), 
Rt<P>  =  (Xt-^)-a1(Xt_1-M)-...-ap(Xt_p-M) 

=  <At(1)-«l>Xt-l  +  (  At<  2  >-a2  )Xt_2  +  ...  +  (  At(  P  >-ap  )Xt_p 

+  Bt-(l-a1...-ap)ji   .  (4.2) 

On  multiplying  (4.2)  by  (Rt_r(P))2  and  taking  expectations,  thus  obtaining 

Cov[Rt(P),(Rt._r(P)  )2]/  we  have  a  sum  of  p  terms  given  by 

P 

Z   E{(At(3)-aj)Xt_j(Rt_r(P))2}  (4.3) 

3=1 
and  a  last  term  involving  Bt  which  is  clearly  zero.   Now  A^  3  )-ou  is 
independent  of  both  the  Xt-j  and  (Rt-r^P^)2,  for  r=l,2,...,  which  may 
nevertheless  be  themselves  dependent.   Thus  the  jth  term  in  (4.3)  becomes 

E(A(i)-aj)E{Xt_j(Rt_r(P))2}  =  0  for  r-1,2,..., 
since  E(A(J))=a-j.  This  completes  the  proof  which  clearly  includes  the 
first  cross  correlation  mentioned  in  the  theorem.  The  proof  highlights  the 
fact  that  it  is  the  independence  of  the  vector  of  coefficients 
{Afc(  *■),  .  .  .  ,At(P),Bt}  on  previous  X-^'s  uhich  creates  the  effect  in  this  type 
of  model. 

The  results  of  the  theorem  can  be  used  to  help  validate  random  coef- 
ficient autoregressive  models;  also  useful  in  this  respect  are  the  non- 
zero higher  order  residual  crosscorrelations,  eg  for  positive  r 
Corr[Rt^P^,(Rt+r^^  )23  •  These  have  been  obtained  for  the  second  order 
autoregressive  exponential  process  studied  in  Tjawrance  and  Lewis  (1984); 
similar  results  for  any  first  order  random  coefficient  model  of  the  type 
(4.1)  are  given  in  the  next  subsection. 

It  is  worth  noting  that  with  G^-r  defined  as  any  reasonable  function  of 
(Xfc_r,Xt_r_i/ . . . ),  a  similar  argument  to  that  given  in  the  proof  shows  that 

Corr[Rt(P),Gt_r]-0,    r=l,2,...  .  (4.4) 
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A  result  of  this  type  is  not,  however,  sufficient  to  establish,  for 
instance,  that  Rt<P>  and  Rt_r<P>  are  independent;  they  are  dependent. 

The  random  coefficient  autoregressive  structure  of  (4.1)  is  not  a 
necessary  condition  for  the  results  of  this  section  to  hold;  this  may  be 
illustrated  by  noting  that  Corr[Rt(P>,(Xt_r-M)2]  is  zero  when 
C12(rVaE{(Xt-M)(Xt_r-ji)2}  satisfy  the  equations 

Ci2(r)=aiC12(r-l)+a2c12<r-2>+-  •  •+«pCi2(*"P>»  r=l,2,...  .  (4.5) 
The  similarity  of  these  equations  to  standard  Yule-Walker  equations  (2.4) 
will  be  apparent.    Thus  any  process  with  this  property  will  have 
Corr[Rt^P^/(Xt-r_R)2]  equal  to  zero  for  non-negative  r. 

4.2  Higher  Order  Dependency  fcr  First  Order  Random  Coefficient 
Autoregressive  Models 

The  model  to  be  considered  is  the  first  order  (p=l)  case  of  (4.1),  now 

to  be  denoted  as 

Xt^AtXt-i+Bt,     t=0,±l,±2, ...  ,  (4.6) 

in  which  At  and  &t  are  independent  within  and  between  each  t;  we  also  write 

a=E(At),   M=E(Xt)=(l-ar1E(Bt). 

Many  basic  mathematical  and  probabilistic  properties  of  this  equation  have 

been  studied  by  Vervaat  ( 1979 ) .   Interest  here  is  restricted  mainly  to  the 

residual  crosscovariances  of  (Rt2,Rt-r)  an<*  [(Xt~/i)2,Rt_r]  where  R^  is  the 

first  order  residual  given  by 

Rt=Xt-^-a(  Xt_  !-/! ) .  (4.7) 

It  has  been  established  in  Section  4.1  that  Cov( Rj- , Rt-r 2  )=Cov(  R^2 , Rt+r  )=° 

for  r=l,2,...,  and  so  now  Cov(Rt2,Rt-r )  for  r=0,l,2,...   is  obtained 

explicitly. 

The  calculation  begins  by  writing  the  required  covariance  as 

Cov(  Rt2 ,  Rt_r  )=e  { [  ( Xt-ji  )-a(  Xt-i-/i )  ] 2  [  ( Xt-r-M  )-a(  Xt^^-n )  ] } 

=C21<  r  )-ac21(  r+1  )-2a{CL11(  r  )-aC1L1(  r+1 ) } 

+a2{C21(r-l)-aC21(r)}  (4.8) 

where,  as  defined  just  before  (4.5), 
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C2i(r)=E{(Xt-ji)2(Xt_r-*0),    C111(r)=E{(Xt-|i)(Xt_1-/i)(Xt_r-^)}. 

The  higher  order  triple  moment  Ci:L1(r)  is  easily  seen  to  be  equivalently 

given  by  aC2i(r-l),  in  the  present  instance,  and  hence  (4.8)  becomes, 

Cov(Rt2,Rt_r)=C2l(r)-ac2i(r+l)-a2{C2i(r-l)-aC2i(r)},   r=l,2,  .  . .   (4.9) 

For  r-0,  with  M3=E[Xt-M-)3]/  there  is  the  seperate  result 

E(  Rt3  )=(  l+2a3  )/x3-3aC21(  1 ) .  (  4 .  10 ) 

The  calculation  of  C2i(r)  is  effected  by  writing  the  defining  equation 

(4.6)  in  the  form 

Xt-M+/x=At(  Xt-x-zx )+( MAt+Bt ) . 
Squaring  both  sides,  multiplying  by  Xt-r-*i,  and  taking  expectations  then 

gives 

C21(  r  )+2/iC(  r  )=a2C21(  r-1  )+2(  /xa2+E(  AtBt )  )C(  r-1 )  (4.11) 

where 

C(r)3Cov(Xt,Xt-r)=ara2,   a2=var(Xt),   a2zE(At2). 

Further  simplifications  of  (4.11)  using  C(r)=ara2  gives  the  recursive 

equation 

C21(r)=a2C21(r-l)+b3ar-1,    r=l,2 (4.12) 

where 

b3=2a2  { jivar(  At  )+Cov(  At , Bt ) } . 

Equation  (4.12)  has  explicit  solution 

C21(r)=a2r^3+b3(a2r-ar)/(a2-a),    r«0,l, ...  .   (4.13) 
Going  back  to  (4.9)  for  Cov(Rt2,Rt-r )  and  using  (4.12),  we  have  finally 
Cov(  Rt2 , Rt_r )=(  a2  +a3a2-aa22-a2 )C21( r-1 ) 

+( l-a2+a3-aa2 )b3ar_1 ,  r=l,2,...  .  (4.14) 
Thus,  (4.14)  and  (4.13)  constitute  the  explicit  solution  for  the 
crosscovariances  of  the  residuals  (Rt2,Rt-r)'  notice  that  when  At  is  the 
constant  a,  the  case  of  the  standard  first  order  linear  model,  these 
crosscovariances  are  correctly  zero.  The  NEAR(1)  model  of  Lawrance  and 
Lewis  (1981)  specializes  (4.6),  by  having  independence  within  each  (At,Bt) 
and  particular  forms  for  At  and  Bt. 

If  the  covariances  (4.14)  are  to  be  converted  to  correlations,  then 
var(Rt)  and  Var(Rtz)  nust  be  obtained.  This  requires  the  following 
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calculation, 

Var(Rt2)  =  E(Rt4)  -  [E(Rt2)]2  (4.15) 

where 

E(Rt2)  =  (1  -  a2)a2  =  Var(Rt),  (4.16) 

E(Rt4)  -  (l+a*)/i4  -  4aC31(l)  +  6a2C22(l)  -  4a3C13(l),         (4.17) 
with 

Ji4=E[  (  Xt-M  )4] ,  Ct  j (  1  )=E[Xt-/i  )i(  Xt_!-M  )3  ]  ,  (  i ,  j  )=(  3  , 1  ) ,  (  3 ,  2  ) ,  (  1 ,  3  )  .  (  4 .  18  ) 

The  joint  moments  in  (4.17)  were  obtained  in  terms  of  the  first  four 

moments  about  the  mean  of  X^,  the  first  three  moments  about  zero  of  At  and 

the  first  two  moments  about  zero  of  B^.   These  are  not  reproduced  here. 

The  explicit  result  for  (4.15)  was  checked  numerically  against  a  slightly 

different   method   of   computations   while   being   used   to   obtain   the 

illustrations  given  in  Figure  1  and  2  for  NEAR(l)  models. 

It  was  remarked  in  Section  3,  following  (3.2)  that  the 
crosscorrelations  of  (X^-M)2  a^nd  Rt+r^*  can  also  be  useful  in  assessing 
higher  order  dependence.  Tn  the  present  case  of  first  order  autoregression 
there  are  the  results, 

Cov[(Xt-ji)2,  Rt_r]  =  C2i(r)  -  a  C21(r+1),  (4.19) 

Var(Rt)  =  (l-a2)a2,  Var[(Xt-ji)2 )]  =  E(Xt-ji)*  -  a*.  (4.20) 

These  formula  can  all  be  applied  to  the  NEAR(  1 )  model  quite  simply, 

noting  that  in  this  case 

E(Atr)  =  o/3r,  E(Btr)  =  rl[p+(l-p)br] 

(4.21) 
b  =  (l-a)0,  p  =  (l-/3)/[l-(l-a)/3]. 

Figures  1  and  2  give  the  results  of  computations  of  Corr[Rt2,Rt-r]  and 

Corr[(Xt;-M-)2/Rt-r]  for  tne  NEAR(  1 )  model.  The  top  left  frame  of  Figure  1  is 

the  linear  KAR(l)  case  for  which  all  cross-correlations  apart  from  lag  zero 

are   zero;   the   other   three   cases   in   Figure   1   each   have   zero 

crosscorrelations  at  negative  lags,  in  agreement  with  the  theoretical 

results  in  Section  4.1,  but  have  some  non-zero  values  at  the  zero  and 

positive  lags.   It  is  evident  that  the  lag  zero  cross  correlations  contain 
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NEAR(1)  RESIDUAL  CROSSCORRELATIONS  (R(T)*2.R(T-R)) 
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Figure  1.  Four  computation  of  the  crosscorrelatlons  for  various  lags, 
between  the  linear  autoregresslve  residual,  Rt^K  and  {Rt(p) }z  for  the 
NEAR(l)  process  ulth  p(  1  ;=a/3  held  constant  at  0.75;  In  effect  the  remaining 
free  parameter  Is  being   varied   through   Its  allowable  range. 


much  of  the  discriminating  information  between  the  four  cases;  this  will  be 
so  for  first  order  autoregresslve  models  in  which  much  of  the  higher  order 
dependency  is  at  lag  one,  and  the  lag  zero  cross  correlation  involves  both 
third  and  fourth  order  aspects  of  the  lag  one  dependency.  For  further 
information  on  the  TEAR<1)  and  PREAR<1)  cases  see  Lawrance  and  Lewis 

(1981). 

Figure  2  gives  Corr[(Xt-j02,Rt_r]  for  the  same  four  cases  used  in 

Figure  1.  The  negative  lags  are  again  zero  for  all  cases,  and  hence  this 
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property  does  not  discriminate  the  linear  from  the  nonlinear  cases.   Rather 
it  is  the  strength  of  the  crosscorrelations  at  positive  lags  which  performs 
this  task,  albeit  less  clearly  than  the  crosscorrelations  of  the  residuals 
and  squared  residuals. 


NEAR(1)  RESIDUAL  CROSSCORRELATIONS  ((X(T)-MU)*2,R(T-R)) 
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Figure  2.  Four  computations  of  the  crosscorrelations,  for  various  lags 
between  the  linear  autoregressive  residual,  R^P) ,  and  (XfU)z  for  the 
NEAR(l)  process  uith  p( 1 )=o@  held  constant  at  0.75;  in  effect  the  remaining 
free  parameter  is  being  varied   through   its  allowable  range. 
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4.3  Higher  Order  Dependency  for  the  Product  Autoregression  Model 

Another  form  of  nonlinear  autoregressive  model,  called  PAR(l),  was 
introduced  by  McKenzie  (1982);  the  basic  idea  of  its  construction  is  to 
consider  the  exponentiation  of  the  standard  AR(l)  equation.  Then  its 
additive  structure  becomes  multiplicative,  and  the  general  form  of  the 
PAR(l)  product  autoregressive  model  model  is 

Xt=(Xt-i)^Bt,    t=0,±l,±2,...,  (4.22) 

where  p  is  the  dependency  parameter,  (Xp<l,  and  {B^}  is  an  independent  and 
identically  distributed  innovation  sequence.  McKenzie  studies  the  model 
when  Bt  is  chosen  so  that  {Xfc}  has  a  gamma  marginal  distribution,  and 
obtains  several  basic  results;  for  instance,  that  p  still  represents  the 
lag  one  autocorrelation  and  that  the  autocorrelations  in  general  satisfy 
the  Yule-Walker  first  order  equations  ((2.4)  with  p=l).  It  is  apparent, 
however,  that  the  linear  conditional  expectation  definition  of 
autoregression  (equation  (2.4))  is  not  satisfied  because  of  the  power  form 
implied  by  (4.22). 

As  with  the  first  order  random  coefficient  models  in  Section  4.2,  we 
consider  the  first  order  residual  Rt,  given  by  (4.7),  and  will  likewise 
determine  Cov(  R^2 , R-t-r )  f°r  r=0,±l,±2, . . . ,  noting  that  for  this  model  these 
correlations  are  non-zero  for  all  lags.  We  will  use  the  general  expression 
(4.8)  in  terms  of  the  third  order  central  moments  C2i(r)  and  C^n(r),  but 
this  time  there  is  no  simple  relation  between  them,  and  both  are  needed  for 
all  lags.  Also,  calculation  of  C2i(r)  and  CA11(r)  must  be  in  terms  of 
their  uncentered  components,  since  these  are  the  quantities  which  can 
immediately  be  determined  from  the  PAR(1)  equation  (4.22).  The  required 
uncentered  moments  will  be  written  as 

e2=E(Xt2),   e3-E(Xt3) 
e21(r)=E(Xt2Xt_r),    e111(r)=E(XtXt_1Xt_r)  (4.23) 

for  r-0,*l,±2, . . . ,  and  there  is  need  to  note  the  special  cases, 
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e2i(0)=e3'   elll(0)se2i(l),   em(l)=e2i(-l)- 
All  these  quantities  will  be  calculated. 

First,  it  is  necessary  to  note  the  following  relations  between  the 
centred  and  uncentred  moments , 

C21(r)=e21(r)-^3-2*iC(r)-M<r2        (all  r) 

Clll(0)=C2l(l),    C111(1)=C21(-1)  (4.24) 

Cin(  r  )=em(  r  )-*i3-^( C(  r  )+C{r-l  )+C(  1 ) } ,   ( r*0, 1 ) . 

The  calculations  of  e2i(r)  and  em(r)  need  to  be  treated  separately 

for  positive  and  negative  lags,  but  follow  in  the  same  general  manner  and 

will  be  illustrated  by  that  for  e2i(r)  for  positive  lags.   By  iterating  the 

PAR<  1 )  equation  (4.22)  r  steps  backward, 

r-1 
Xt  =  Xt_rP(r)  n  Bt_iP(i> 
i=0 

where  p(r )spr,r=0,l, . . . .   Squaring  this  equation,  multiplying  it  by  Xt_r 

and  taking  expectations,  gives 

r-1 
e2i(r)=E{Xt_r2P(r>+1}  nEJBt-i2^1)}.  (4.25) 

i=0 

To  obtain  the  expectations  in  the  repeated  product,  taking  the  2p(i)th 

power  of  (4.22),  leads  to 

E{Xt2p<i>}=E{Xt_12PP(i)}E{Bt2P(i)}, 

and  hence  dropping  the  unnecessary  suffix  t,    to  the  result 

E{B2p(i)}=E{x2p(i)}/E{x2p(i+l)}> 

Now  (4.25)  can  be  expressed  purely  in  moments  of  X,  as 

e2i<r)=E{X2P(r>+1}E(X2)/E{X2P(r>}.  (4.26) 

To  proceed  further,  invoke  a  gamma  marginal  distribution  for  X  with  density 

f(x)=e^3_1e-ex/r(/3),   e,/3>o,  x>o;  (4.27) 

this  has  mean  /3/e,  variance  /3/e2,  third  central  moment  2/3/e3,  and  there  is 

the  kth  moment  result  E(  Xk  )=r(  |3+k  )/©kr(  /3 ) .    From  (4.26)  and  similar 
expressions  we  then  have 
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*    ,<r-\  f/3(/3+l)(0+2pr)/e3,  r=l,2 

62l(r)    ~    U0fpl*l>«HPl*l+l>/e3.  r=-l,-2,...  (4-28) 

6lll(r)  "   U/9+plr')(/94-p+plrl+l)/e3#        r=-l,-2,...  (4>29> 

with  e21(0),    eni(O)   and  e1^11(l)  being  given  by  the  special  cases  of 

(4.23).      Use  of  (4.24)  gives   finally  the   required  expressions 

[2/3pr/e3  r=0,l 

C2l(r)  "   Up»*l+p2lrl)/e3,  r=-l,-2,...  (4'30) 

r         /r\  f/3(pr+pr+1)/e3,  r=2,3 ,AO,* 

Cm(r)  -  l/3{(1+3)p|r|+l+p2|r|+l_flp|r|-l}/e3   r=1/2,...  <4'31> 

All  these  may  be  used  in  (4.8)  to  obtain  the  desired  Cov(Rt2,Rt_r )  for 
r=0,±l,±2, . . . ;  also  the  simpler  covariance  or  correlation  of  (X^-M)2  and 
Rt-r  is  easily  available  from  (4.19)  and  (4.31). 

Conversion  of  these  covariances  to  correlations  requires  Var(Rt2)  as  at 
( 4 . 15  ) ;  the  required  intermediate  results  are 

M4  =  9  33/©*, 
C31(l)  =  30(0+2)p,  C13(l)  =  /3{3/3+(l+p)(2+p)}p, 

C22(l)    -   02+2/3p+2/3(  /3+2  )p2  .  (4.32) 

These  then  give  the  explicit  expression 

Var(Rt2)    =    {(  9/3-1  )/32   -    4/3(/3+6)p2    +    12/3p3 

+  /3(9/32-/3+16)p*-12/3ps-4/3p6}/e*  (4.33) 

Together  with  Var(Rt)  which  is  simply  (l-p2)/3/e2,  Cov(Rt2,Rt-r )  can  ^ 
converted  into  Corr(Rt2/Rt-r )• 

Figures  3  and  4  illustrate  the  residual  and  squared  residual  cross 
correlation  for  the  PAR( 1 )  model.  Figure  3  shows  how  these  residual  cross 
correlations  vary  over  the  range  of  p  values  from  p=0  to  p=0.9  for  an 
exponential  marginal  distribution.  At  p=0,  the  PAR(l)  is  an  TDD  process 
and  the  residuals  are  trivially  independent,  resulting  in  the  only  non-zero 
crosscorrelation  at  lag  zero.  As  the  p  value  increases  the  dependence 
spreads  out,  most  strongly  at  lag  minus  one  and  positive  lags.  An 
interesting  feature  is  the  lag  zero  crosscorrelation  which  changes  from 
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being  strongly  positive  at  p=0  to  approximately  zero  at  p=0.75,  to 
moderately  negative  at  p=0. 

PAR(1)  RESIDUAL  CROSSCORRELATiONS  (R(T)*2.R(T-R)) 
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Figure  3.  Computations  of  the  residual  crosscorrelations ,  Corr(Rf,z ,R^.r) 
for  the  PAR(l)  model.  Since  /3=i  the  marginal  distribution  is  exponential. 
The  lag  one  serial  correlation  is  increased  from  p(l)=0.0  (upper  left)  to 
p(l)=0.9   (loner  right). 


Figure  4  gives  four  different  gamma  cases  of  the  PAR( 1 )  residual  cross 
correlations,  all  with  p=0.75.  The  gamma  shape  parameter  0  takes  the 
values  0.5,  l.o,  2.5  and  10.0;  these  cases  indicate  that  changes  in  the 
gamma  parameter  cause  only  modest  changes  in  detail  of  the  cross 
correlations. 
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PAR(1)  RESIDUAL  CROSSCORRELATIONS  (R(T)*2.R(T-R)) 
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Figure  4.  Computations  of  the  residual  crosscorrelatlons ,  Corr(R-^z ,R-^.r) 
for  the  PAR(l)  model.  The  p  for  each  case  of  0.75,  and  the  figures 
Illustrate  the  effect  on  the  residual  crosscorrelatlons  of  changing  the 
Index  (3  of   the  gamma  distribution   through   0.5,    1.0,    2.5  and   10.0. 


5 .   CONCLUSIONS 

A  methodology  for  analyzing  higher  order  dependence  in  nonlinear  time 
series  with  pth  order  autoregressive  correlation  structure  has  been 
proposed.  It  utilizes  standard  uncorrelated  linear  autoregressive 
residuals,  and  the  crosscorrelation  function  of  these  residuals  and  their 
squares.  The  behaviour  of  this  crosscorrelation  function  has  been  utilized 
for  two  rather  different  types  of  nonlinear  model:  random  coefficient 
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autoregression  and  multiplicative  autoregression:  the  behaviour  has  been 
shown  to  allow  discrimination  between  models  in  the  same  class  with  the 
same  marginal  and  autocorrelation  structures. 

The  residuals  crosscorrelation  function  provides  a  partial  analysis  of 
third  order  information  in  the  time  series;  it  does  not  attempt  to  capture 
all  third  order  information,  which  is  the  aim  of  such  techniques  as 
bispectral  analysis  and  which  will  often  be  intractable  with  nonlinear 
models.  Being  based  on  standard  linear  residuals,  the  analysis  extends 
rather  than  replaces  conventional  residual  analysis. 

Developments  of  the  analysis  which  focus  on  the  directionality  implicit 
in  many  time  series  are  being  investigated;  reversed  residuals  assume  a 
reversed  directionality  and  allow  exploration  of  the  consequences  of  such 
an  ass sumption. 
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